Exponential Decay of Correlations in a Model for Strongly Disordered 2D Nematic 

Elastomers 
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Lattice Monte-Carlo simulations were performed to study the equilibrium ordering in a two- 
dimensional nematic system with quenched random disorder. When the disordering field, which 
competes against the aligning effect of the Frank elasticity, is sufRciently strong, the long-range cor- 



relation of the director orientation is found to decay as a simple exponential, e 



-r/e 



The correlation 



length 5 itself also decays exponentially with increasing strength of the disordering field. This result 
represents a new type of behavior, distinct from the Gaussian and power-law decays predicted by 
some theories. 
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The effects of quenched random disorder are crucial for 
physical systems as diverse as random ferromagnets, vor- 
tex lattices in type-II superconductors, structural glasses, 
foams and liquid crystals. It is therefore of great impor- 
tance to attain an understanding of the way in which the 
range of the ordering of such systems varies as a function 
of the relative strengths of the ordering and disordering 
fields. In this Letter we present results from a study of 
one simple model system. The behavior we observe is 
qualitatively different from any that has previously been 
discussed. 

Strong disorder, manifesting itself in stochastic values 
of certain microscopic parameters, does not permit the 
system to achieve its ordered ground state at low tem- 
peratures. This results in a variety of glass-like struc- 
tures with only short-range order, but with correlations 
distinctly different from those in a high-temperature liq- 
uid state. An example is the Sherrington-Kirkpatrick 
spin glass 0J^, which is represented theoretically by an 
Ising model with random exchange constants, the vari- 
ance of whose fluctuations acts as a measure of the de- 
gree of disorder. When the characteristic energy of the 
disordering field is small, or the number density of dis- 
crete random impurities is low, the system is capable of 
establishing ideal ordering on a small but macroscopic 
length scale. On larger scales the cooperative effects of 
the random field may prevail over the demand for uni- 
form thermodynamic order, and the long-range correla- 
tions decrease with distance in systems with sufficiently 
low space dimensionality. The possibility then exists of 
separating the effects of elasticity, which on a continuum 
level penalizes the deviations from the ground state or- 
der, from those of the continuum, coarse-grained, weak 
random field. 

The concept of a well-established local order that is 
distorted on length scales larger than a characteristic co- 
herence length to form textures with zero global av- 
erage requires a multi-component order parameter. The 
simplest example is the vector order parameter of fer- 



romagnets, in which a dilute random array of impuri- 
ties gives rise to correlated spin-glass textures ||,^. A 
system that is in many ways analogous to the ferromag- 
net is a nematic liquid crystal having a uniaxial tensor 
order parameter, where the principal axis (the director 
n) can form birefringent textures. Nematic elastomers 
are very close in their behavior to the correlated spin 
glasses, with network crosslinks acting as local sources 
of weak quenched disorder |^ . A recent analysis based 
on the concept of hierarchical replica symmetry breaking 
Q has sugested that at long distances a quasi long-range 
order in vortex arrays pertains, i.e. the fluctuations di- 
verge logarithmically and a simple power-law decay of 
long range correlations is predicted. This is in contrast 
to the general arguments of Imry and Ma who pre- 
dict that the long-range correlation of a continuous order 
parameter should decay as exp [— (f'/^d)'' ''j with d the 
dimensionality of space. 

All the existing continuum models are based on the 
notion of weak disorder, which is difficult to control para- 
metrically. The region of applicability of such a concept 
could be tested experimentally, if a series of similar sys- 
tems with systematically increasing random field were 
available. In reality, however, random magnets exist only 
with a limited range of possible impurities [Q, and it is 
in general difficult to monitor the spatial variation of lo- 
cal magnetization in textures. Nematic elastomers offer 
more promise for experimental studies because the bire- 
fringent textures are readily accessible to several stan- 
dard optical techniques |^ . This gives us the motivation 
to study such systems in depth. 

In this Letter we use the direct lattice Monte-Carlo 
simulation of equilibrium textures in a model correspond- 
ing to nematic elastomers, which thus shows strong sim- 
ilarity to spin glasses. We look for the change in long- 
range correlations that occurs in response to a variation 
of random field strength. Because we need to analyze 
large systems, for computational efficiency we consider 
only the case of = 2. The nematic director field n(r) is 
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assumed to be of constant (unit) amplitude and its direc- 
tion to vary in the plane so that a single angle 9(x, y) is 
the only coordinate to enter the Hamiltonian. This cor- 
responds to a locally well-established orientational order 
forming a large-scale texture, and is closely related to a 
2D XY model {d = 2,n ~ 2), the difference lying in the 
invariance of the nematic elastomer to inversion of the 
director n to — n. 

We consider a nematic with misoriented rod-like 
crosslinks, which locally impose a preferred orientation 
direction on the director field. This is analogous to a 
continuum ferromagnet with randomly distributed impu- 
rities. The Hamiltonian of such a system can be written 
as: 



(1) 



Here K is the Frank elastic constant (in the single- 
constant approximation) and 7 the strength of the ran- 
domly oriented field, whose local direction is defined by 
the unit vector at the crosslink location £. In a typical 
elastomer the crosslinks are separated by distances much 
smaller than the expected correlation length S^di and typ- 
ically of order 1-5 nm. 

To perform the Monte Carlo study we first make the 
assumption that the essential element of the disorder lies 
in the orientational randomness of the crosslinks rather 
than their spatial distribution. This allows us to use a 
model in which the sites i form a regular lattice, each 
site of which contains a random field k(^), with no corre- 
lations between neighbors. The density of crosslinks, po, 
is thus assumed uniform. 

The minimization of the free energy (|l|) determines in 
each cell the local director orientation n{re), which at 
zero temperature will be a function only of the dimen- 
sionless parameter D* = Aj/K. The resulting discrete 
dimensionless Hamiltonian, scaled by the Frank constant 
K, takes the form 
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where S are the nearest-neighbor lattice vectors, and the 

— 1/2 

lattice spacing, pg , becomes the unit of length. 

The key information we seek is the form of the cor- 
relation function C(r) that describes the distance over 
which the director remains largely unchanged in orien- 
tation. In the systems we shall study, in which the di- 
rector is confined to a 2-dimcnsional manifold, this re- 
quires computation of the average over positions ro of 
C(r) = 2(cos'^[0(r + ro) — 9{ro)]) — 1. For a macroscopi- 
cally isotropic (polydomain) system we expect this quan- 
tity to be a function only of |r|, and so we have chosen 
to evaluate it as the average over all directions of r. 

Rather than perform the simulations on a very large 
system we studied lattices of modest size (typically 60 x 
60) and performed many repeated evaluations of C(r) 



for different realizations of the random array of crosslink 
orientations. In this way it was possible to estimate the 
reliability of the results from the width of the spread in 
C(r) produced. The effect of the lattice size was found 
by comparing the results for lattices of 60 x 60 and 40 x 40 
sites. The values of C(r) were indistinguishable for small 
r|, and differed by only a few per cent even when |r| was 
half the size of the smaller lattice. 

The procedure adopted to find the minimum energy 
configuration was to run the Metropolis algorithm for the 
Monte Carlo method at a finite temperature and then to 
reduce the temperature to zero for the remainder of the 
run. This process of annealing and quenching avoided 
having the system be stranded in some of the smaller lo- 
cal energy minima. Some further tests of the robustness 
of the results were made by varying the initial conditions. 
A starting point of complete alignment, in which the di- 
rectors at each lattice site were initially mutually parallel, 
produced results for C(r) that tended to be larger than 
when the starting point was a completely random set 
of orientations. Each series of runs was then continued 
until the results from the two opposing starting points 
converged onto each other. 

A more physically significant choice to be made con- 
cerned the boundary conditions used. The evaluation of a 
reliable form for C(r) is greatly facilitated by using peri- 
odic boundary conditions, since then all lattice points are 
equivalent. On the other hand, the absence of any phys- 
ical boundary to the system has the consequence that 
topological defects cannot evaporate at a surface, and can 
only be eliminated by a process of mutual annihilation. 
The weakness of the effective attraction that occurs be- 
tween certain types of defects hinders this process, and 
makes some defects very long lived. The presence of a 
free boundary thus greatly speeds the process of equili- 
bration, as defects are able simply to wander off the edge 
of the system. The disadvantage that comes with the 
use of free boundaries, however, is the fact that not all 
lattice sites are equivalent, and the calculation of cor- 
relation functions must be limited to pairs of sites that 
are in the interior of the system. In the simulations to 
be described below, it was indeed found that the use of 
free boundary conditions led to a greatly accelerated ap- 
proach to equilibrium, and so this scheme was generally 
employed. 

The decision to adopt free boundary conditions was 
reinforced by an examination of the consequences of ap- 
plying a boundary condition intended to represent the 
effects of surface anchoring. In this modification, a uni- 
form field was applied only to the lattice sites lying at 
the boundaries, and the simulation was run with the 
two opposing initial configurations, the totally random 
and the completely aligned. It was then found that the 
two sets of correlation functions never approached each 
other, no matter how long the run was continued, the 
initially aligned state giving a larger value for the corre- 
lation function at all separations r than the one starting 
from totally random director field. To avoid this effect 
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we have preferred the free boundary conditions. 

The correlation function C(r) was calculated as the 
average of 2(n£p • n^(,+^)^ — 1 over all sites I contained 
within a circular shell of radius r, centered on site 
of thickness one lattice spacing. The calculations thus 
describe the effects of the competition between the ran- 
dom field, which tends to destroy all correlations, and the 
nematic interaction, which tends to introduce long-range 
correlations in the director orientation. 
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FIG.l The orientation correlation function C(r) is shown as a 
function of distance, measured in lattice spacings, for values of 
D* , the relative strength of the random field due to crosslinks, 
equal to 0.5 (diamonds), 1.0 (crosses), and 3.0 (squares). 

When D* is large, the effect of the random fields of the 
crosslinks dominates, and the decay of C(r) is extremely 
rapid. As D* is reduced, the range of the correlation 
function gradually increases. In the range studied, for 
which 0.5 < D* < 3.5, the correlation function is found 
to approximate fairly closely a simple exponential. This 
is illustrated in Fig. , which shows on a semilogarithmic 
plot the form of C(r) for D* = 0.5, 1.0, and 3.0 with 
lines representing the exponential fit. 

The necessity of checking the robustness of the calcu- 
lation by choosing two opposite starting configurations 
was confirmed by performing an average of 50 runs of 75 
000 time steps at each of the 3600 sites, starting from 
both ordered and random initial states. The two results 
were in disagreement by about 30% at large r. Only af- 
ter the runs had been continued to 150 000 time steps 
did the two curves approach each other sufhciently to be 
indistinguishable. 

The linearity of the plots in Fig. tells us that the cor- 
relation function can be written in the form C(r) = e~^^^ 
over the range of D* studied. This form for the correla- 
tion function permits us to identify ^(D*) as the correla- 
tion length for this system, measured in units of the lat- 
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tice spacing^Q ' . We expect it to be related to the "do- 
main size" |3). When we examine the variation of S,{D*) 
with D* we find that this in turn is fitted by an expo- 
nential to within the precision of our computations, as is 
shown in Fig. , and thus is of the form £,{D*) oc . 
This result is somewhat surprising in view of the fact 
that ^ must tend to infinity as D* vanishes, and suggests 
that there are two distinct regimes for the correlation be- 
havior in this randomly disordered nematic system. It is 



unfortunate that the amount of computation needed to 
produce reliable results increases very steeply as D* is 
reduced below 0.5, and that it was consequently not fea- 
sible to perform simulations to follow the form of £,{D*) 
down into this small-D* crossover region. At the other 
extreme, where D* becomes larger than 3, the correlation 
length diminishes to a fraction of a lattice spacing. (This 
statement is not as devoid of meaning as it appears, since 
it is merely a way of saying that the correlation in orien- 
tation between nearest neighbors has fallen to a number 
less than e~^.) 
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FIG. 2 The effective correlation length ^ is found to vary expo- 
nentially with the strength of disorder D* . 

Further calculations were made for a non-zero applied 
external field, which enters the Hamiltonian as — (h-n^)^ 
and introduces an additional dimensionless parameter 
{h'^/K) in Eq.(||). The results here were less unexpected, 
in that the correlation function C(r), after an initial rapid 
decay, approached asymptotically a constant value that 
reflects the susceptibility of the system and was, accord- 
ingly, proportional to {h? /KY for small h. 

To summarize, the results of our simulation describe 
the equilibrium ordering of nematic elastomers, whose 
large-scale alignment is disturbed by the quenched mis- 
orienting effect of network crosslinks. The continuum 
coarse-grained model used to describe this system is 
closely similar to that for a spin glass. We considered the 
case of a two-dimensional XY model, which is a good de- 
scription of a thin flat sample of nematic elastomer with 
the director confined to this plane by planar boundary 
conditions. This is the typical experimental arrangement 
described earlier. Accordingly, the model falls into the 
class of (c? = 2, n = 2) continuum Hamiltonians with the 
quenched random field weak enough to allow formulation 
in terms of a continuum displacement field 9{x,y). 

The classical arguments of Imry and Ma 0] predict 
correlations to decay as C(r) ^ exp[— (r/^2)^i while a 
treatment accounting for replica symmetry breaking |^ 
predicts C{r) ^ {flC)~^i with x a positive fractional 
exponent. Our results appear to be quite unambigu- 
ous and to provide a totally different correlation law: 
C(r) w exp(— r/^) with the characteristic decay length 
^ ~ 23exp(-1.16L'*) whenever D* > 0.5. This behavior 
cannot persist down to vanishingly weak disorder, as it 
fails to predict an infinite value for ^ when D* 0. This 
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suggests that there is a crossover to another regime at 
smaller values of D* = j/K than considered here, and 
which might be within the range of validity of the ana- 
lytical predictions |^,||. The fact that an extrapolation 
of our data to D* ~ gives a decay length roughly equal 
in size to the linear dimensions of the lattice used in the 
calculation is almost certainly only a coincidence; the 
large- 1?* behavior from which the extrapolation is made 
should be independent of the lattice size. 

We are left with a number of unanswered questions. 
What is the physical explanation of the exponential 
form of the decay length? Over what range of disor- 
der strength does it persist? Is this result unique to two 
dimensions or does it have a more general applicability? 
The resolution of these problems will require some deeper 
insight and a significant dedication of computational re- 
sources. 
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